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ABSTRACT 



Context. High-resolution numerical methods have been developed for nonlinear, discontinuous problems as they appear in simulations 

of astrophysical objects. One of the strategies applied is the concept of artificial viscosity. 

Aims. Grid-based numerical simulations ideally utilize problem-oriented grids in order to minimize the necessary number of cells at 

a given (desired) spatial resolution. We want to propose a modified tensor of artificial viscosity which is employable for generally 

comoving, curvilinear grids. 

Methods. We study a diff'erential geometrically consistent artificial viscosity analytically and visualize a comparison of our result to 

previous implementations by applying it to a simple self- similar velocity field. We give a general introduction to artificial viscosity 

first and motivate its application in numerical analysis. Then we present how a tensor of artificial viscosity has to be designed when 

going beyond common static Eulerian or Lagrangian comoving rectangular grids. 

Results. We find that in comoving, curvilinear coordinates the isotropic (pressure) part of the tensor of artificial viscosity has to be 

modified metrically in order for it to fulfill all its desired properties. 

Key words. Hydrodynamics - Radiative transfer - IVIethods: numerical 



> 

On 
<N 

O 

o 



^ 



1. Introduction 

In astrophysics a multitude of systems and configurations are 
described with concepts from hydrodynamics, often combined 
with gravitation, radiation and/or magnetism. Mathematically 
radiation hydrodynamics (RHD) and magnetohydrodynamics 
(MHD) are described by systems of coupled nonlinear partial 
diflTerential equations. The Euler equations of hydrodynamics, 
the Maxwell equations as well as radiative transport equations 
are hyperbolic PDEs that connect certain densities and fluxes via 
conservation laws. The numerical solutions of these equations 
essentially need to comprise this quality. Today there exists a 
wide range of numerical methods for conservation laws that en- 
sure the conservation of mass, momentum, energy etc. if applied 
properly, and multiple fields in physics and astrophysics have 
adopted these sophisticated numerical methods for studying var- 
ious applications. 

Standard numerical schemes for partial diflTerential equations 
are established under the assumption of classical diflferentiabil- 
ity. Routine finite difference schemes of first order usually smear 
or smoothen the solution in the vicinity of discontinuities as they 
come with intrinsic numerical viscosity. Standard second order 
methods show something to the eflTect of the Gibbs phenomenon, 
where oscillations around shocks emerge. In the past decades so 
called high-resolution methods have been developed in order to 
achieve proper accuracy and resolution for nonlinear, discontin- 
uous problems as they appear also in RHD or MHD. One of these 
strategies is the concept of artificial viscosity which we will also 
briefly motivate in sub section |1.1| 



In higher-dimensional problems this artificial viscosity 
emerge s as a tensorial quantity. We will demonstrate this in sub- 
The result we want to present in this paper can be 



1.3 



section 

seen as a tensor analytical consequence of the artificial viscosity 
in general curvilinear coordinates when using consistent metric 
tensors. In section [2] we will propose a correction for the com- 
monly used tensor of artificial viscosity for curvilinear grids. 

This correction is motivated by astrophysical applications 
where one considers comoving nonlinear coordinates repre- 
sented by non-conformal (non-angle preserving) maps from 
spherical coordinates. The authors are currently investigating the 
generation of grids that are asymptotically spherical but which 
allow certain asymmetries that can be found in rotating configu- 
rations, nonlinear pulsation processes etc. This new approach to 
grid-based astrophysical simulation techniques will be addressed 
extensively with numerical applications in a future paper. 

As an example of non-conformal two dimensional coordi- 
nates, Figure[T]shows a grid that corresponds to the map (x, y) -^ 
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a^TT^ - I6a2^ + a2a37r^^ 
47r(l+a2^) 



a^rj^j sin rj 



1 + 



(1) 



which gives the polar coordinates for the choice of parameters 
(ai,a2,<23) = (1,0, 0). In such a nonorthogonal grid the metric 
tensor is no longer diagonal and one has to consider a consistent 
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(d) 



Fig. 1. Example of a non-conformal non-steady 2D grid with oblate- 
ness governed by three choices of parameters (ai,a2,a2) as defined in 
equation ([T]). 



differential geometric approach to the formulation of the govern- 
ing equations of RHD and MHD, and also to the mathematical 
formulation of the artificial viscosity, which will be stressed in 
Section 121 

The benefit of the consistent formulation can be seen when 
we consider time-dependent grids, e.g. when using time- 
dependent parameters (ai,a2,a3) in ([T]). We refer to the Ap- 
pendix[5]for the depiction of the system of equations of RHD for 
generally comoving curvilinear coordinates with time-dependent 
metrics correspondingly. 

1.1. Brief Introduction to Conservation Laws 

For the sake of stringency we recapitulate some important results 
from the theory and numerics of conservation laws which can be 
found e.g. in (LeVeque 1991 ) or ( [Richtmyer and Morton|1994] ). 

The equations of RHD and MHD form a system of hyper- 
bolic conservation laws that describe the interaction of a density 
function d (x, : ^"^ X [0,oo) -^ W and its flux f (d) : W -^ 
j^mxn Equation ^ shows how a concrete choice for the density 
and the flux field can look like in a given coordinate system. 

The temporal change of the integrated density in a connected 
set Q c R" then equals the flux over the boundary ^Q, i.e. 



dt fd dV -F ff • n d^ = for aU t > 0, 



(2) 



where n is the outward oriented normal of the surface. 

The system is called hyperbolic if the Jacobian matrix V^f 
associated with the fluxes has real eigenvalues and if there exists 
a complete set of eigenvectors. In case of MHD and RHD this 
property has a direct physical relevance ( [Pons et al.|20Q0| ). 

Assuming f to be a continuously diflTerentiable function, 
equation ^ can be rewritten via the divergence theorem as 



//<^"' 



-h divx f(d)) dV d^ = for aU ^ > 0, D c R^ (3) 



t Q. 
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which gives the system of partial diflferential equations for the 
density function d: 



dtd -h divx f (d) = for all r > 0, x g : 



(4) 



With an initial condition d(x, 0) = do(x), x g R", this is called 
the Cauchy problem. 

In order to illustrate the connection of hydrodynamical appli- 
cations to this formalism, we express the Euler equations in the 
form (|4]). The appearing variables are the gaseous density p(x, t), 
the gas velocity u(x, t), the inner energy 6(x, t) and the gaseous 
pressure tensor P(x, t). Considering the diflTerential form (|4]), we 
recognize the continuity equation, the equation of motion and the 
energy equation as the components of the hyperbolic problem. 
In case of the most relevant problem, that of 3D hydrodynamics, 
the density and its flux are given as 



(p) 
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f(d) = 


puvJ + P 
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[pen + (Puf) 



j,5x3 



(5) 



Eor a given coordinate system with base vectors e/, the ten- 
sorial fields are given explicitly as (using the Einstein notation) 



f 



d = 






f(cl) = 



pw e^ 






(6) 



The gaseous pressure tensor can be assumed to be isotropic in 

most applications, which means that P = g^^p where p(x, t) is 

.J . 
the scalar gas pressure and g^-^ = e^ e^ the contravariant metric 

tensor. In case of adaptive grids respectively moving coordinate 
lines, the base vectors are time-dependent as well, i.e. then C/ = 
e/(x, t). 

Since even the simplest examples of one-dimensional scalar 
conservation laws like the Burgers' equation have classical solu- 
tions only in some special cases, one has to broaden the consid- 
ered function space of possible solutions. For the so-called weak 
solutions, we appeal to generalized functions where the discon- 
tinuities are defined properly. The generalized concept of diflTer- 
entiation of distributions shifts the operations to test functions 
7 : R" X R^ D G ^ R (G open) which are infinitely diflTeren- 
tiable and have a compact support (meaning that for each y there 
exists a closed and bounded subset K such that y(x, t) = for 
all X e G\ K). We denote this space of test functions by D(G). 
In this generalized space of solutions the Cauchy problem ^ is 
written as 

r r (dtd -h divx f(u))r dV d^ = for afl y g D(G). 

Jt>o Jw 

The weak formulation of the conservation law ^ is obtained by 
shifting the derivatives to the test functions by partial integration, 
and by using the compactness of the support. We get that the 
following has to hold for each y e D{G): 



// 



{Ad,y + m)^^j)AVdt 



/ 



r(x,o)do(x)dy. (7) 



t>0 



The function d g L'^ is called a weak solution of the PDE (|4]), 
if it satisfies ^ and d ^ U with do g L!^ . However, there is 
a small drawback. This weak solution is not necessarily unique 
and usually further constraints have to be imposed in order to 
guarantee its uniqueness. This leads us to the actual topic of this 
paper. 
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1.2. Introduction to Artificial Viscosity 

For most physical problems it is naturally sufficient to look for 
weak solutions from the function space of piecewise continu- 
ously differentiable functions. Constraining the space of solu- 
tions in this way, we call the physical variables d weak solutions 
of the Cauchy problem (|4]), if they are classical solutions wher- 
ever they are continuously differentiable, and if at discontinuities 
(shocks) they satisfy additional conditions in order to be physi- 
cally reasonable (we elaborate on these conditions below). 

The mathematical theory provides several techniques to dis- 
tinguish physically valuable solutions out of a manifold of math- 
ematically possible. One method is to add an artificial viscosity 
term to the right hand side of (|4]), to get the equation: 



(9,d-hdivxf(d) = £yAd, s>0 



(8) 



and then consider the limiting case s ^ 0. This idea is motivated 
by physical diffusion which broadens sincere discontinuities to 
differentiable steep gradients at the (microscopic) length scale of 
the mean free path of the particles. The physical solution of the 
weakly formulated problem is thus the zero diffusion limit of the 
diffusive problem. However, in practice this limit is difficult to 
calculate analytically, and hence simpler conditions have to be 
found. A common technique to do this is motivated by contin- 
uum physics as well. Here an additional conservation law is set 
to hold for another quantity - the entropy of the fluid flow - as 
long as the solution remains smooth. Moreover, it is known that 
along admissible shocks this physical variable never decreases, 
and the conservation law for the entropy can be formulated as an 
inequality. 

We denote the (scalar valued) entropy function by cr(d) and 
the entropy flux function by 0(d), and they satisfy 



dtcr(d) -h divx 0(d) = 0. 



(9) 



Assuming the functions to be differentiable, we may rewrite this 
conservation law via the chain rule and the equation ^ as 



Vdcr(d)divxf(d) = divx0(d), 



(10) 



where in higher-dimensional case the appearing matrices of gra- 
dients have to fulfillfurther constraints, see e.g. (G odlewski and| 
Raviart|1992| ). For scalar equations, it is always possible to find 
an entropy function of this kind. Furthermore it is assumed that 
the entropy function is convex, i.e. 



V^o- > 0, 



for all deU. 



(11) 



To get our actual entropy condition, we first rewrite our entropi- 
cal conservation law ^ in the viscous form 



dtcr(d) -h diVx 0(d) = 6:VdCr(d)Ad. 



(12) 



Integrating over an arbitrary time interval [^o, ^i] and a connected 
set ri c R", and using partial integration, we find that 



(dtcr(d) -h diVx 0(d)) dV dt 



ff 

to Q 

h 

=£ f r(VxdVdO-(d))-nd^ d^ 

to do. 

ti 

-.//v.,dVV5d)V„<,ciVd. 



(13) 



to Q. 



>0 



When we now consider our non-diffusive limit s ^ 0, the first 
term on the right-hand side vanishes without further restriction 
whereas the second term has to remain nonpositive. With partial 
integration and divergence theorem we get our entropy condition 



fcr(d(x,ti))dV< fcr(d(x,to))dV 



t\ 



// 

to da 



(14) 



0(d) • n d^ dt. 



For bounded, continuous pointwise solutions d* of (12 ) such that 
d* -^ dfov s ^ 0, the vanishing viscosity solution d is a weak 
solution of the initial value problem ([3]) and fulfills entropy con- 
dition ( p4| ). Generally spoken, applying the entropy condition to 
systems with shock solutions unveils those propagation veloci- 
ties that ensure that no characteristics rise from discontinuities 
which would be non-physical. For detailed motivation, strin- 
gent argumentation and proofs to mathematical techniques pre- 
se nted in this sectio n we refer to ( |Hartenetal.|1976| ) respectively 
to ( |LeVeque]T99T] ). 



1.3. Numerical Artificial Viscosity 

As mentioned we are looking for high-resolution methods for 
nonlinear PDEs derived from hyperbolic conservation laws. In 
the past decades major eflforts have been made in developing nu- 
merical methods for these problems that are at least of second or- 
der. One patent attempt to finding such a high-resolution method 
is to adapt a well-known high-order method for linear problems 
for nonlinear problems (such as the Lax-WendroflT scheme ([Lax 
[aSdW endrofl^ I960)). 

As illustrated above we can add an artificial viscosity term to 
the conservation law in a way that the entropy condition is sat- 
isfied and non-physical solutions are excluded. We are keen to 
design this viscosity in such a manner that it aflfects sincere dis- 
continuities but vanishes sufficiently elsewhere so that the order 
of accuracy can be maintained in those regimes where the solu- 
tion is smooth. The idea of numerical artificial viscosity was in- 
spired by physical dissipation mechanisms and dates back more 
than half a century to (von N eumann and Richtmyer|195Q) . 

We denote as customary the approximate solution of the ex- 
act density d(x, t) at discrete grid points d(xj, tn) by D" and set 

D = [Di . . . Dy^]^, where k is the total number of grid points. 
The numerical representation of the flux function f (d) is denoted 
respectively by F(D), where [F(D)]y = f(Dy). The numerical 
flux function gets modified by a an artificial viscosity Q [D]^ for 
instance in the following way: 

[Fvisc(D)],- = [F(D)],. - h [Q(D)] . (D,-^i - D,). (15) 

where h denotes the size of the spatial discretization. Since the 
original design of this additional viscous pr essure in the scalar 



form Q = C2p(Au)^, C2 e R as suggested in ( von Neumann and 



|RichtmyeT|| 1 950 ) for one dimensional advection dfd -\- adxd = 
Qdxx^, it has undergone a number of modifications and gener- 
alizations. It has tu rned out t o be numerically preferable to add 
a linear term (see (Landshoff 1955)) in order to control oscilla- 
tions. Generalizations to multi-dimensional flows mostly retain 
the original analogy to physical dissipation and reformulate the 
velocity term accordingly, see e.g. (Wilkin s 1980 ). 

The artificial viscosity broadens shocks to steep gradients at 
some characteristic length scale, but should not cause too large 
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smearing. The concrete composition and implementation of this 
artificial viscosity coefficient Q depends on the application. As 
an example we discuss the following form of the tensor of the ar- 
tificial viscosity in higher-dimensional RHD numerics. Similar 
forms of artificial viscosity can be found also in pure hydrody- 
namics and MHD calculations in 2D and 3D. 

Tscharnuter and Winkler (Tscharnu ter and Winkler] |1979| ) 
have pointed out that the viscous pressure in 3D radiation hydro- 
dynamics has to unravel normal stress, quantified by the diver- 
gence of the velocity field and shear stress, which is expressed 
by the symmetrized gradient of the velocity field according to 
the general theory of viscosity. It is designed to switch on only 
in case of compression (divx u < 0), and this is all ensured by 
the form 

Q = -^24cPmax(-divxU,0)HVu]^ - -edivxuj, (16) 

where the symmetrization rule is defined componentwise for the 
lower indices as 



([Vu],X-^-= -(ViUj-^VjUi). 



2. Numerical Artificial Viscosity in Curvilinear 
Coordinates 

We introduced artificial viscosity in form of a three dimen- 
sional viscous pressure tensor ([16]) in section [O] along the lines 
of ("Tsch arnuter and Winkler||1979| ). In this section we want to 
point out, how such a definition must be adapted for curvilinear 
coordinates in order to ensure tensor analytical consistency. 

When formulating PDEs derived from hyperbolic conserva- 
tion laws on a curvilinear grid, the tensorial equations ^ have to 
be transformed to the according coordinate system. Not only the 
vectorial and tensorial quantities have to be transformed but also 
the differentiation operators, in particular the divergence opera- 
tor in our case. The appropriate framework to do this is provided 
by differential geometry. Like the gradient of a scalar is natively 
a covector, there are rules for co- and contravariant indices of 
tensoLS such as the one we are interested in. The crucial term 
in ( 16 ) is the symmetrized velocity gradient Vu that accounts 
for snear stresses, and one sees that the form ( [16| ) comes into 
conflict with the demand of vanishing trace (TrQ = 0) when the 
divergence term is simply of the form e divx u , as we find it 
commonly in several MHD and RHD grid codes. 

Proposition: The correct form of the viscous pressure tensor 
( p^ in general coordinates is 



Q = -4tscP max(-divx u, 0) [Vu] 



^gdivxU 



(17) 



We show that this tensor has the desired properties. The viscous 
pressure tensor must be symmetric by definition, i.e. Qij = Qjt, 
which can be easily verified from (Vf) . Also, the trace of theten- 
sor has to vanish (TrQ = Q\ = o\ as pointed out in |von Neu- 
mann and Richtmyer|1950| ). To show that this holds for (pT|), we 
first consider its native covariant components 



Qij = -/i2 max(-divx u, 0) I -(V/t/y -\- V jUt) - -gtj divx u | , (18) 



where we have renamed q\l^^^^p = J^i- Next, we need to raise an 
index with the metric, 

&j = Qijg' 

= -ii2 max(-divx u, 0) -g^X^iUj -\- V jUi) - -d'j diVx u| , 

and use the essential identity g^^gij = g^j = S^j. The Ricci 

Lemma Vtgjk = digjk - T^tjgik - ^^ikgji = for the fundamental 
tensor naturally also holds for the contravariant components and 
we can permute g into the derivatives Vig^^Uj and V jg^^ui which 
yields twice the divergence V/w^ = diVx u when we conduct the 
contraction j -^ i. In three dimensions the summation S^i = 3 
and we obtain our desired result 



1 



-(2 diVx u) - -3 diVx u I = 



The commonly used ( see e.g. (|Dorfifl 999) in RHD, ( [Iwakami 
et al.|20Q8| in MHD, ( jFryxeU et al^p OQQ) in MHD) form of Q 
(|16|) is not compatible with these requirements since the sym- 
metrization is only defined for lower indices, whereas the unit 
tensor e of a metric space is only defined for mixed indices, 
meaning there is no such thing as Sjj. However, the above men- 
tioned and other authors such as ( [Mihalas and Mihalas]|1984| ) 
have neglected that little inconsistency since they have con- 
sidered mixed indices from the start respectively cartesian or 
affine coordinates. Nonlinear corrections have been suggested 
by ( [Benson and Schoenfeld||199"3] ) albeit they do not explicitly 
concern curvilinear coordinates and are based on a TVD ap- 
proach. 

In several hydro- and MHD-codes that include non-Cartesian 
grids such as Pluto ( [Mignone et al.|2Q(J7] ), the geometric source 
terms are coded explicitly for several geometries (polar, cylindri- 
cal, spherical), and not only for the artificial viscosity flux. The 
suggestions made e.g. by (Marcel and Vinokurp 974 ) lead to ge- 
ometrical source terms that correct curvilinear grid effects. How- 
ever the strong conservation form as elaborated by ( Warsi|1981| 
would need to appeal to our differential geometrically consistent 
approach in order to deal with the viscosity in an intrinsically 
consistent way. Especially when the metric tensor itself is not 
only a function of space but also time-dependent (as discussed in 
Section [T]), the latter approach reaches its limits. Our correction 
affects curvilinear coordinates in multiple dimensions, whereas 
it is not necessary that the coordinates are orthogonal respec- 
tively the metric tensor does not need to be diagonal. Our initial 
motivation to study more general coordinates comes from the 
idea to generate problem-oriented coordinate systems for astro- 
physical numerical calculations. In a following paper we want to 
present some feasible approaches to grid generation under cer- 
tain physical restrictions. Such nonlinear grids that are adaptive 
in multiple dimensions have time-dependent metric tensors and 
thus benefit directly from our consistent definition. On the con- 
trast, when using adaptive mesh refinement, the metric tensor 
remains geometrically constant in time. 

In order to support the theoretical results in this work, in the 
upcoming section we study as an example a very simple velocity 
field with non- vanishing divergence and visualize the according 
artificial viscosities for the two presented cases. 



3. Application and Visualization 

The most common application of curvilinear coordinates in 3D 

is the map (x,y,z) ^ (r g R+, ?^ g [0,7r], (p e [0,2;:]) with 
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X = rsin?^cos^, y = rsin?^sin^ and z = rcos§ as spher- 
ical coordinates. The corresponding diagonal covariant met- 
ric components in this simple orthogonal case are given by 
diag(l,r, rsin?^). 



3.1. Toy Model Velocity Field 

As the presented considerations for artificial viscosity on non- 
steady curvilinear coordinates originate from astrophysical ap- 
plications, we want to consider a velocity field with a certain 
practice in RHD. We study a toy model of a self- similar fluid 
flow solution, namely the velocity field given by 



UEx 



^Jx^ -\-y^ -\- z^ ^ 



(19) 



here in Cartesian coordinates. Such self- similar solutions ap- 
pear in idealized spherical models of stars for example as shocks 
driven by radial stellar pulsations. 

This vector field is obviously symmetric with respect to 
the origin and has a non- vanishing divergence divx Uex = 2/r. 
The covariant components of this vector field are given in any 
other coordinates by scalar product with the base vectors, i.e. 
^Ex,/ = Uex • e/. This leads to the covariant components (1,0,0) 
in spherical coordinates. The nonzero covariant components of 
the tensorial part ([Vu]^ - |edivxU) of the artificial viscosity 
( p^ are given for this field by 



Qrr = -^, 

3r 



Qoi 



r sin 6 



(20) 



In the following section we visualize the tensor of artificial vis- 
cosity (TAV) for the velocity field ([19]). A uniform distribution 
of the leading eigenvalues over the whole domain is expected 
due to the symmetry of the vector field. 

One easily verifies the identity Q\ = summing over the 
mixed components Q\ = -2/3r, Q% = l/3r, Q^^ = l/3r. 



With the previous version of the TAV from ( Tscharnuter and 
Winkler 1919) we would get the following covariant components 
of the artificial viscosity tensor. 



Qrr = - 



3r' 



Qee = 



2 
3r 



= -hrsin^^. (21) 

3r 



The visualization of this non-metric version of the artificial vis- 
cosity for the symmetric velocity field ([19]) shows obviously 
a field unequal in strength and direction over the whole do- 
main. In a numerical calculation this will clearly lead to ar- 
tificial anisotropics in the flux of the density field and destroy 
all eflTorts in constructing a higher-order conservative numerical 
scheme with artificial viscosity. 

3.2. Visualization of scalar and tensor fields 

We can see the incorrect vs. correct behavior immediately by 
even just displaying the major eigenvalue or trace as one indi- 
cator. However, since the major eigenvalue represents only one 
degree of freedom out of the six available in the tensor field, a 
technique depicting all six components is a more objective way 
for validation. 

We used the Vish Visualization Shell ( [Benger et al.|2007| to 
numerically sample d2Q| and ( [2T] ) on a uniform grid for analyz- 
ing scalar fields (Fig.[2JFig.[3]) and a radial sampling distribution 
(Fig. |4]) for the full tensor field. 

Fig. [2] displays a structure of the eigenvalue corresponding 
to the major eigenvector of the TAV on the XZ plane (i.e., in 



i 




(a)i. 



I 




(b)^. 



Fig. 2. Major eigenvalue the incorrect (left) and correct (right row) 
viscosity, shown along the XZ plane. 



the plane y = 0), evidently showing some asymmetric 'funny' 
coordinate-dependent behavior of the incorrect TAV, whereas the 
correct TAV is radially symmetric, as desirable. The tensor field 
is symmetric and of rank two, however it is not positive definite 
and may exhibit vanishing trace. The correct TAV is trace-free in 
the entire domain, whereas the trace of the incorrect TAV ranges 
through positive values (for large radial distances) to large neg- 
ative values close to the coordinate origin, as depicted in Fig. [3]. 
A direct visualization method depicting the full six-dimensional 
degrees of freedom the tensor field is thus favorable, or even 
rather essential. However, many direct visualization methods 
for tensor fields require positive definiteness and are thus not 
applicable to this data. Only very few methods are suitable for 
general tensors. Fig. |4] shows so-called Reynold glyphs ( [Moore 
[eraL] ri995) for the TAV field. A Reynold glyph is the surface 
generated by mapping a tangential vector v(P) at each data sam- 
ple point P as 



P^Q(v(P),v(P)). 

Such glyphs shown at each sampling point provide a direct vi- 
sualization of the full six-dimensional degrees of freedom of 
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I 




(a) TrQi, 




(b) TrQeorr 

Fig. 3. Trace the incorrect (left) and correct (right row) viscosity, 
shown along the XZ plane. Note the range of values in Fig.|3b]- what we 
see is just numerical noise, the trace is perfectly zero within numerical 
precision. 



the tensor properties. Reynold glyphs are able to also de- 
pict negative definite tensors, whereas a quadric surface (ellip- 
soids representing P ^ 1/ ^/Q (v(P), v(P)) becomes hyperbolic 
for negative eigenvalues and problematic for visualization pur- 
poses. The Reynold glyph directly show the "directional value" 
Q (v(P), v(P)) of the tensor field Q in the direction v around a 
the sampling point P - the resulting surface is intersecting the 
sampling point P whenever Q(v(P),v(P)) = 0, which is the 
case for points where the tensor is degenerating and not posi- 
tive definite. In such areas the glyph will visually appear like 
two intersecting surfaces corresponding to the isopotential sur- 
faces of second order spherical harmonic functions. These both 
surface components represent positive and negative eigenvalues 
of the tensor field - if both positive and negative component 
"counter-balance" themselves they therefore indicate vanishing 
trace, which is the sum of the eigenvalues. If only one surface 
component is visible, then the tensor is either positive or negative 
definite on that certain point. 

We used a radial sampling for the direct visualization of the 
tensor field in order to minimize coordinate artificats. As de- 



picted in Fig.|4b]and Fig.[5b]for the correct TAV all "modes" are 
equivalently represented, indicating vanishing trace of the tensor 
field while being radially aligned with the underlying coordinate 
system. In contrast, the incorrect TAV, Fig. |4a| and Fig. |5a| ex- 
hibits a dominantly negative trace. 
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Fig. 4. Reynold glyphs ('Moore etal.'1995\ of the incorrect and correct 
viscosity tensor. The incorrect tensor showing a strongly negative com- 
ponent, whereas the correct tensor is balanced, indicating zero trace. 
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(a) Detail of the incorrect viscosity 




(b) Detail of the correct viscosity 



Fig. 5. Detailed Reynold glyphs ( Mo ore et al.|l995| ) of the incorrect 
and correct viscosity tensor. Glyphs of the correct tensor show spherical 
harmonics that are balanced in their positive and negative half- widths. 



4. Conclusions 

We studied a generalization of the tensor of numerical artificial 
viscosity for curvilinear coordinates and compared our result to 
previous definitions found in literature. We analyzed a sym- 
metric toy velocity field and visualized its viscosities. Clearly, 
the non-metric version of the TAV as used by many authors of 
hydro- (HD) respectively magneto-hydro- (MHD) or radiation 
hydrodynamic-codes (RHD) leads to incorrect results in curvi- 



linear coordinates, whereas our suggestion for the numerical ar- 
tificial viscosity gives geometrically consistent results. 

5. Appendix 

As mentioned in section [T] the benefit of the strong formulation 
of Cauchy problems arising from physical conservation laws will 
be discussed in more detail here. 

Following the ideas of ( [Warsi|198T] ), in non- steady coordi- 
nates the geometrically consistent strong form of tensorial con- 
servation laws ([2]) is given as 



dt{^\d)^di{^\m)'e^) = 0, 



(22) 



where d and f are to be decomposed according to their native ten- 
sorial components. The scalar multiplication with the /-th con- 
travariant base vectors e^ yields a projection on the contravariant 
coordinate lines which in case of generalized grids can differ in 
direction and length of their covariant counterparts. Equation 
( |22| ) gives also the integral form of the conservation law, which 
should be treated numerically in a correct way for non-steady co- 
ordinates in any finite- volume discretization. The important dif- 
ference to the componentwise structure V(.) = d(.) + r(.), where 
ChristoflTel symbols account for the geometry is that in this case 
undifferentiated terms arise (see ( [Marcel and Vinokur||1974| )) 
which act like geometric sources in the equations and destroy 
conservativeness. A comprehensive proof of Vinokurs theorem 
using differential forms can be found in ( [Bridges |2QQ8| ). 

With non- steady curvilinear grids not only the nonlinearity 
of the metric tensor but also its time-dependence has to be taken 
into account numerically. The motion of the grid itself and its 
implications on the formulation of the set of equations respec- 
tively the calculation of the occurring fluxes is discussed in the 
following section. 

5.1. Adaptive Grids 

In fluid dynamics we distinguish two main reference systems 
that suit unequally for various applications. The Eulerian frame 
is the fixed reference system of an external observer in which the 
fluid moves with velocity u whereas the Lagrangian approach 
describes the physics in the rest frame of the fluid. Between 
these two systems, the transformation of an advection term for a 
density d (that moves with a relative velocity u) is given via the 
material derivative Dtd = dfd + u • Vd. 

Hence, when we work with comoving frames, the coordinate 
system respectively the computational grid is time-dependent. 
There is a number of purposes where strict Eulerian or La- 
grangian grids are suboptimal and thus we need to consider the 
generalized the concept of the comoving frames. 

We want to evolve the strong conservation form for time de- 
pendent general coordinate systems. The time derivative of a 
density d in the coordinate system S(^) relative to a (e.g. static) 
coordinate system S(q,) is given by ^?d(^) = ^?d(Q,) + V(Q,)d ^?X(^) 
and from the view point of system S(q,), the time derivative is 
given as 



dtd(^) = dtd, 



iP) 



■ V(,)dT, 



(23) 



where x denotes the grid velocity. The second term on the right 
side we call grid advection, see e.g. (Warsi"1981) and (Thomp-[ 
sonet al.,1985j . An inhomogeneous advective term of a conser- 
vation law K = d^ + div (u^d) (in a fixed coordinate system) is 
given in the case of moving grid as 



K = d, - X 



Vd^ + div(u^d). 



(24) 
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When we apply the appropriate transformation prescriptions to 
the spatial derivatives, we gain the following form 



^—di ( ^e^d) + -^di ( ^u . e^d) . (25) 



K = dt-x 



Now this is not yet geometrically conservative, since it is not of 
an integral structure. Following the idea of the Reynolds trans- 
port theorem we consider the temporal derivative of the volume 
respectively the determinant of the time dependent metric tensor 
^J\g(xJ)\ in order to study the conservation of a density function 
in variable volumes and obtain the strong conservation form for 
time dependent coordinate systems: 

^K = d, ( ^d) + dt ( ^e^" . (u - x)d) (26) 



For the full analytic derivation we refer to ( [Thompson et al.| 
[1985 ) again. Defining the contravariant velocity components rel- 
ative to the moving grid by U^ = e^ - (u- x) the above equations 
yield 

^K = d, ( ^d) + dt ( ^^^d) . (27) 



5.2. Set of RHD Equations in Strong Conservation Form 

We exhibit the system of equations of radiation hydrodynam- 
ics in somewhat simplified formulation. The following system 
has been basis for a number of implicit RHD computations (see 
e.g. (Dorfi^J_999)). All the astrophysical assumptions, implica- 
tions and simplifications can be found in ( [Mihalas and MihaJas 
1984|). In this paper we only want to emphasize the structural 



form of such a set of equations in a strong conservation form for 
comoving curvilinear coordinates. Note that for scalar equations 
the only eff'ectively remaining geometric term inside the deriva- 
tives is the volume element Ylgl- The vectorial equations contain 
however also the time dependent base vectors. 

5.2.1. Continuity Equation 

The strong conservation form of the continuity equation 

dtP + div (pu) = 
is given for time depended coordinates by 

5.2.2. Equation of IVIotion 

The equation of motion that we consider in radiation hydrody- 
namics contains the conservation of moment of plain fluid dy- 
namics ([5]), the radiative flux as a coupling term (H, with the 
specific Rosseland opacity of the fluid kr), gravitational force 
(G) and the artificial viscosity term (Q from ([17])): 

dtipu) -F div (pu^u -h P -h Q) -h pG KRpU = 0. 

c 

We investigate the elements of the energy stress tensor a little 
closer before we give the consistent strong conservation form. 
We define an eff'ective tensor of gaseous momentum R that ac- 
counts for the motion of the coordinates as 

R = r^^e/Cy = p(u - x) u^. 
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The isotropic gas pressure tensor P is defined by the scalar pres- 
sure and the metric tensor as P = pg. The viscous pressure ten- 
sor Q is to be modified in the way we suggested in ^Tl} . Since in 
most applications of RHD with self-gravity involved, the grav- 
itational force G is determined by solving the Poisson equation 
for the potential O, we substitute G = -VO. The equation of 
motion in strong conservation form is then written as 



^r(A^pu) + ^.(A^(R + P + Q)-e^) 



47T 



^pdi ( A^Oe^) - —KR ^pH = 0. 



The ^th component of the strong conservation equation of mo- 
tion is given by the ^th Cartesian component of the unit vector, 
e.g. in spherical coordinates e^ = cos if sin ?^e^ -h sin ip sin ?^e^ -h 
cos ?^e^ and its derivatives. The projection of each physical ten- 
sor on the contravariant coordinate lines can be simplified by 
its contravariant components with respect to its covariant basis 
without losing strong conservation form, i.e. f • e^ = f'^^j. We 
prefer this form with contravariant components since it meets 
the native design of the stress tensor R and the pressure tensor 
P. The tensor of artificial viscosity as given in ([18]) is brought to 
contravariant form by the summation 



Qlmg g 



li fyinj 



li mj 



;1 



^y(V/w^ + V,-wJ-g^^-divu 



and then the equation of motion as 



dt ( ^|\g\pu^et) ^^t{^f\^\ {r^^ + p^^ + &^) ej) + 



47T 



pdi ( ^Oe') - —KR ^\pW^i = 0. 



(28) 



5.2.3. Equation of Internal Energy 
The energy equation 

dtipe) -F div (up6) -F P • Vu^ - 47iKpp(J - S) -\-Q'Vu^ = 

accounts for the thermodynamics of the fluid, namely the energy 
balance including kinetic and pressure parts as well as inner en- 
ergy. Latter is a thermodynamic quantity which is associated 
with the equation of state. The specific inner energy (e) is in 
case of an ideal fluid its thermic energy. Another term comes 
from the energy exchange with the radiation field ((/ - S )-term) 
containing the specific Planck opacity kr and viscous energy dis- 
sipation, expressed by the contraction of the viscosity with the 
velocity gradient Q • Vu^. 

Since we assume isotropic gas pressure F = pgwe can re- 
formulate its contribution via the Ricci Lemma and obtain a very 
simple scalar expression. 

P . Vu'^ = g'^pViUj = pViu' = /7divu 



The viscous energy dissipation is given by the contraction 
Q . Vu^ = a^ViUj = a^ {diUj - Y\jUk) =: Ec 



-diss? 



and the strong conservative form of the energy equation is then 
given by 

^t{^|\g\pe)^^i{^\pee^ • (u - x)) + ^J\^\pdiYU- 
- An ^\Kpp{J -S)+ ^Ed,ss = 0. 
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5.2.4. Equation of Radiation Energy 

We write a simplified frequency integrated radiation energy 
equation in the comoving frame as follows 

dtJ + div(u/) + cdivH + K • Vu^ + cxp{J -S) = 0. 

For the scalar energy input of radiative pressure into the material 
we define a new coupling variable 

and in strong conservation form the equation of radiation energy 
is given by 

dt ( A^/) + di {^W- (/(u - X) + cH)) + 

+ A^Pcoup + MCXP{J - ^) = 0. 

5.2.5. Radiative Flux Equation 

Another vectorial conservation law, the radiative flux equation, 
remains to be written: 

dtH + div(uH) + cdivK + H . Vu'^ + cxr^ = 0. 

We define an eflTective radiative flux tensor L analogously to the 
eflTective tensor of gaseous momentum: 

(u-x)H=:L = f^e/ey 

and for the contribution of radiative momentum to the material 
H • Vu^ we define another coupling variable F with components 

The geometrically conservative form of the radiative flux equa- 
tion in non- steady coordinates is then written as 



dt{^\li) + dt{^\e^-(L + cK)) + 
+ ViglFcoup + ^|\g\KRpli = 0. 



coup 
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